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When a steady temperature gradient is established be- 
tween the two ends of a piece of solid bar, heat current 
will flow from high to low temperature end. According 
to Fourier's law of heat conduction the current density 
is proportional to the temperature gradient and mathe- 
matically it reads as 



J(x) = -kVT(x), 



(1) 



where the constant of proportionality k is known as the 
thermal conductivity of the solid. Conduction of heat in 
solid by its very nature is a non-equilibrium process. This 
is an area of Physics, where the idea of non-equilibrium 
statistical mechanics can be applied to in order to find 
the underlying physical conditions for the validity of this 
law in case of solid. Various numerical and analytical 
studies confirm that the heat transport in one dimen- 
sional system exhibits anomalous [l| behaviour. It means 
that the thermal conductivity for such a system is not 
found to be an intrinsic property of the material. It 
shows a power law dependence k ~ N a , where N be 
the linear size of the system. There are studies on differ- 
ent models which predict divergent (0 < a < 1) thermal 
conductivity [lHH . There are also some oscillator models 
that give non-divergent (a < 0) thermal conductivity 0] 
in one dimension. The anomalous behaviour of thermal 
conductivity is also observed in two dimensional system. 
Numerical study indicates a logarithmic divergence of 
thermal conductivity k ~ In N. A power law behaviour Q 
is also observed in such a system. 

There are strong numerical evidences [9( that indicate 
the validity of Fourier's law of heat conduction in one and 
two dimensional systems with pining and anharmonicity. 
An extensive investigation on heat transport in a three 
dimensional disordered harmonic crystal has been carried 
out recently [lOj. The numerical simulation indicates the 
normal transport of heat when this system is subjected to 
an external pining potential. Though it is not been ver- 
ified numerically, but a finite conductivity is predicted 
for this disordered system from analytical arguments. A 
more recent simulation study 11| establishes for the first 
time the validity of this law in three dimensional anhar- 
monic crystal. It thus also establishes the fact that the 



process of heat conduction in three dimensional geometry 
is diffusive in nature. Apart from bringing in a tempera- 
ture dependent contribution to the thermal conductivity, 
which is indeed the case for real systems, it is confirmed 
that anharmonicity provides a condition which is suffi- 
cient for normal heat transport in a solid. In this letter 
we give an exact analytical derivation of Fourier's law of 
heat conduction in three dimensional harmonic crystal. 
We find that in the continuum limit the thermal conduc- 
tivity is finite and does not depend on the system size. 

We consider a cubic crystal in three dimension. The 
form of the Hamiltonian 



H 



(2) 



The displacement field x n is defined on each lattice site 
n = (rii,n 2 ,n 3 ) where ni = 1, • • • , N, n-z = 1, • • • , W 2 , 
and 113 = 1 , ■ • • , W3 . Here e denotes the unit vector 
in the three directions. We choose the value of mass at- 
tached to each lattice point and the harmonic spring con- 
stant as one. We have Langevin's type heat baths that 
are coupled to the surfaces at n\ = 1 and n\ = N and 
are maintained at temperatures Tl and Tr {Tl > Tr) 
respectively. Hence the equation of motion of a particle 
at the site n reads 

X n = ~ XX Tn ~~ X n+e) - l(S ni .l + S nuN )x n 
e 

+ (<W^ + <W^). (3) 

We have chosen the noises to be white and they are un- 
corrected at different sites. Noise strength is specified 
by 



(vZ R (t)v% R (t)) = 2lT L , R 6(t - 0<Sn,n', 



(4) 



where we have chosen the Boltzmann constant fcg = 1. 
We use the periodic boundary conditions for the displace- 
ment field and the noises in n 2 and n 3 directions: 



x a +(o,w 2 ,o) (t) 

^n+(0,VK2,0)(*) 



X n (t) — I n +(0,0,W 3 ) (*) 



(5) 
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These periodic boundary conditions lead to the following 
expansion of x a (t) and r]^' R (t): 



x a (t) 



¥W H E »m (P2,fls, t)e l(p2n2+P3n3)a , 



VW2W3 



(6) 

a 2(p 2 n2+p3«3)a 

} 

(7) 



where a be the lattice constant of the crystal. Upon 
substitution of Eqn.© and (J7J into Eqn.© we obtain 



Vj = -VjkVk - jWjkVk + fj 



(8) 



where 



W jk = Sj,iSk,i + Sj,Nfa,N, (9) 

fj(P2,P3,t) = Sj t lf L (p2,P3,t) + 5j, N fR(P2,P3,t)(10) 

the N x N matrix 

(2ul -1 ... \ 



V 



-1 2w 2 -1 
-1 2ul -1 



V 



(ii) 



-1 2w 2 / 



and 



w 2 ( P2 , P3 ) = l + 2sin 2 (^)+2sin 2 (^). 



Here j, k — 1, • • • , AT. We have also assumed here that 
ya(P2,P3,t) = = y N+ i(p 2 ,P3,t). To solve Eqn.© we 
diagonalize the matrix V. The solution of the N order 
equation \ V — a 2 l\ = gives the eigenvalues of V as 



a 2 k (p 1 ,p 2 ) = 2w 2 (pi,_p 2 ) + 2cos 



kn 
N + l 



(13) 



The j-th component of the normalized eigenvector corre- 
sponding to the eigenvalue a 2 is given by 



(-l) J+1 sm' ' 



at + i 



AT+ 1 



(14) 



The diagonalizing matrix A thus reads as Ajk — Oj such 
that A T A = / and j4 t W1 = a 2 , where (a 2 )^ = a 2 <^- fe . 
We introduce a new set of coordinates £j as 



% (P2 , P3, *) = (Pa , P3, *)• 



(15) 



The equation of motion of £j in matrix form can be writ- 
ten as 



(16) 



where the symmetric matrix Z = A T WA, and / = A T f. 
In the steady state limit (t >> 1/7) we are interested in 
the particular solution of the set of equations of motion 
of £. We use the Fourier transform of 



/oo 1 />oo 7 

— e,(a;)e-* and = / f£/;Me* 
-00 27r J 2ir 



(17) 



in Eqn. (|16l) and obtain 



(-w 2 <5j fe + (XjSjk + iju>Z jk )€k(u) = fj{w)- (18) 

Since the dynamics of the system in the steady state is 
governed by the noises, we decompose £j (w) as 



(19) 



and then using this decomposition into Eqn. (jT5]) we ob- 
tain 



b(L0) 



lu 2 — a 2 — ijuj 



(20) 



Now upon substitution of Eqn. (fT9|) into (fT7f along with 
the use of Eqn.([l0]), ([20]) we obtain 



t;j{P2,P3,t) = - — 
J -00 2tt 



up- — a 2 — ijoj 



x[a^ } fL(p2,P3,^) + a ( N /fl(P2,P3,w) 



0) 



(21) 



(12) Now the use of Eqn.flT}, dTUJ) and (JTTJ) into © gives 



(fL,R{P2,P3,u)f L M(p 2 ,p 3l UJ )) 

= 4%jT LiR S(u + u}') S P2+p > 2t0 5 P3+p > 3t0 . (22) 

To compute the correlation between position and veloc- 
ity we use Eqn. (|2"Tj) and (|2"2"|) and after performing a fre- 
quency integration using delta function obtain 



(£,k 1 (P2,P3,t)^k 2 (P2,P3^')) 

2 7 (a[ kl) a[ k2) T L + a^a^T R ) I c (t - t') 



(23) 



where 



Ic(t-t') 



duj 



j(t-t') 



2ir (w 2 — a 2 ^ — i^u))(uj 2 — a 2 ^ + i^uj) 
Performing the integration over uj we obtain 
e — kI*— 



Ic(t - t') 



4A d (/3 1 ,/3 2 ) 

+i<(t-f )*(*'-*)] 



[I>(t-t')6(t-t') 



(25) 
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where 

I>(t-t' 



) = (cos^i - cos/3 2 ) 2 



-7 2 (2cjo + cos/3i + cos/3 2 ), 
2(cos/?i — cos /3 2 ) cos(wfc 1 |t — t'\ 



(26) 



-{(4^o 2 



■ 3cos/3i + cos/3 2 ) 



_ _7 

xsin(w fcl |i-f'|)}, 

/3i - cos/3 2 ) cos(wfc 2 |t - i'Q 



(27) 



I<{t-t') = 2(008,., ..... ...... . . 

- — { {4uj 2 + cos Pi + 3 cos /? 2 ) 

xsin(a; fc2 |<-<'|)}, 

7T*1,2/(JV+1), 



01,2 



*1,2 



- 774. 



(28) 
(29) 

(30) 



It is clear that I c (t 
when t = t' 



t') -> 0, when |t — i'| — !• 00 and 



_ COS Pi - COS ffg 

MUj ~ 2A d (ft,/3 2 ) ' 



(31) 



For 1 < \ki ~ k 2 \ < iV — 1, / c (0) remains finite when 
N tends to infinity. According to Eqn. (fT4")) the fac- 
tor appeared in Eqn.([23]l (a[ kl) a[ k2) T L + cl^Tr) = 
2 (T L + ( - 1 ) k 1 +k2 T R ) sin /3i sin p 2 / (N + 1 ) . It implies that 
even for zero momentum modes (p 2j 3 = 0), which appear 
owing to the periodic boudary conditions imposed on the 
displacement field in n 2 and directions, the equal time 
correlation in Eqn. (|2"3")l goes as N~ a (1 < a < 3) when 
TV —> 00. The fall of this correlation as a negative power 
of N in the thermodynamic limit indicates that the bal- 
lastic transport remains absent from the conduction pro- 
cess of heat[8j. 

Heat current density j n from the lattice site n to 
n + e 1; where e ± = (1, 0, 0), is given by[l[ 



(32) 



The average heat current density per bond|ll| 



N-l W 2 W 3 



J = 



2W 2 W S (N- 1) 



ijEEE J*- 



(33) 



ni— 1 712 — 1 713 — 1 



We substitute Eqn.© and (fT5j) in J and after performing 
the summations over n 2 and 713 obtain the average heat 
current density per bond in the steady state limit as 



1 



N N-l 



2W 2 W 3 (N 



it E E E (°£ 



P2.P3 fel,fe2 = l ni = l 

x (tffij + (pa.Ps, (-P2, -P3, *)>• 



(34) 



We now use Eqn. (fT4|) to evaluate the sum 

N-l 

\ a n 1 + l "m )\ u n 1 +l' u ni 

ni = l 

= 2(1 - (-l) fel+fc2 )sin/3isin/3 2 
1 



1 



and then using 



J = — 



cos p 2 — cos Pi 
23]) and (HH) obtain 

2 1 (T L -T R ) 
(N +1) 2 (N -1)W 2 W 3 



(35) 



A' 



E E 

P2,P3 fci,fc 2 =l 



x(l - (-l) fc i+ fe s) 



fe.+fc^ sin 2 ^ sin 2 /3 2 
Ad(A,/3 2 ) ' 



(36) 



The factor (1 — (— l) fe i+ fc 2) ensures that the summation 
over ki and /c 2 will be non zero only when ki+k 2 is an odd 
number and hence we take the factor (Tl + (— l) kl+k2 TR) 
out of the summation as (Tl — Tr). In the continuum 
limit, when a — > and W 2 . 3 — > 00 keeping a W 2i3 at fixed 
values, we convert the discrete sums over p 2 and p 3 into 
integrals: 



E 



aW 2 . 
2tt 



dp2,: 



(37) 



Evaluation of the integrals [12j over p 2 and p 3 gives 

2 7 (T L -T fl ) 



J = 



where 



I(N,j) 



1 



JV — 1 



JV 



(38) 



(iV + 1) 2 



]T (i-(-i) fei+fe2 ) 



Here the function 

A(Pi,p 2 ) = (cos/?i -cos/3 2 ) 2 

+7 2 (6 + cos/?i +cos/3 2 ). (40) 

I(N, 7) is zero if fci and fc 2 simultaneously take even in- 
teger values or odd integer values. Assuming that N be 
an even number and using the fact that the summand of 
Eqn. (|39[) is symmetric in respect of the interchange of Pi 
and P2, we rewrite the double sum of 



I(N, 7) = 



(7V+1) 2 



N / 2 -23-25 

sin Pi sin p 2 

A(A,/3 2 ) 



/ (i 1 l;(4 7 2 /A(/3 1 , / 3 2 )) 2 j, (41) 
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where (3i = 2-Kji/(N + 1) and (3 2 = 7r(2j 2 - 1)/(N + 1). 
Again in the continuum limit we convert this double sum 
into integrals. In this limit a — > and N — > oo keeping 
Na at a fixed value. Defining the integration variables 
in this limit as 6*1,2 = 27r/i j2 /(7V + 1), we convert the 
discrete sums into integrals: 



N ■ 



N/2 

31,2=1 



-> - / d9 h2 



(42) 



I(N, 7) thus takes the form 



3(7) = lim I(N,j) 

TV— >oo 







d6 1 / d9 2 



sin Vi sin t'2 







F[\, 1, 1;(4 7 2 /A(0 



A(0i,0 2 ) 
)) 2 



(43) 



Hence we obtain the steady state current density per 
bond in the continuum limit 



J = — n 



(T L - T R ) 
N - 1 ' 

where the conductivity 

K = 273(7). 



(44) 



(45) 



Here k is found to be independent of the size of the sys- 
tem. The variation of the thermal conductivity 
function of 7, as given by Eqn. (|45|) . is plotted in Fig|T] 
Here 7 appears as a constant in the dissipative force 



0.2 0.4 



FIG. 1: (Color online) Plot of k as a function of 7 

term of the Langevin's equation. Physically this force 
term denotes a viscous force experienced by the particles 
of Brownian like at the boundary surfaces of the crys- 
tal owing to collisions with the particles of fluid which 



seems to constitute the heat baths |l3|. The increase of 
7, reduces the mobilities of the Brownian particles and 
thereby reducing their velocities TH 14 [ . Consequently, 
the velocities of the particles at the surfaces next to the 
boundaries will also fall because those are connected by 
springs with the particles at the boundaries. This fall of 
velocities of the particles at the neighbouring surfaces of 
the boundaries will reduce the rate of flow of heat from 
the boundaries to the crystal itself and thereby reducing 
the thermal conductivity of the system. Hence, it jus- 
tifies reasonably the nature of variation of k with 7 as 
shown in FigfTJ 

The average of the square of velocity of a layer at n\ 
reads 



^avg 



W 2 W 3 

n 2 = l "3 = 1 
N 

W ■•]]';■ ^ 53 



1 



1 



E 



"ill Til 



P2,P3 fcl,fe 2 = l 

x (P2,P3, t)ik 2 (-P2, -Pa, t)). (46) 

We use Eqn. (f2"Tj) to compute the velocity- velocity corre- 
lation as 

(61 (P2,P3, t)ik 2 (-P2, -Pa, t)) 
= -^j(T L + (-l) k ^)sm(3 lS m(3 2 
2luq + cos j3i + cos /?2 



(47) 



Upon substitution of Eqn. (|47j) into Eqn.(|46| and evalu- 
ation of p 2 and P3 sum in the continuum limit along n 2 
and na directions, give 



v 2 av Jn l ) = h L (n u N)T L + h R (ni,N)T R (48) 



where 



h L (n u N) = 



h R (n u N) 



N 

E 



(TV +1)2 A. A(/3i,/3 2 ) 

— A 



(49) 



A(/?i,/3 2 



x sin(ni/3i) sin(ni/? 2 ) sin/?i sin/3 2 , 

x sin(ni /3i ) sin(ni /3 2 ) sin f3i sin f3 2 , (50) 
{(cos ,Si - cos/3 2 ) 2 

x[l-F(l/2,l/2,l;(4 7 2 /A(/3 1 ,/3 2 )) 2 )]} 
+7 2 (6 + cos/3! +cos/3 2 ). (51) 



Our evaluation suggests that for 7 = 0.01, Hl tends to 
0.0396 and and Iir tends to and 0.0396 at n\ = 1 and 
m = N respectively when N — > 00. It indicates that as 
Hl and hji are monotonically decreasing and increasing 
functions of ri\ respectively, i; 2 attains a minimum at 
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N = 500 



100 200 300 400 

n. 



FIG. 2: (Color online) Plot of «^„ 9 as a function of m 



N = 700 



100 200 300 400 500 600 700 800 

FIG. 3: (Color online) Plot of v 2 avg as a function of ni 



any layer in the region between m = 1 and n\ = N and it 
is also evident from our plots given in Figj2]and[3] Since, 
v avg( n i) i s proportional to T(ni), the temperature of the 



layer at n\, T(ni) also exhibits a minimum in the region 
1 < n\ < N. This concave upward nature of T(?ii) has 
also been predicted in Ref. jFU] 

In summary, we have given an exact analytical deriva- 
tion of Fourier's law of heat conduction in a three di- 
mensional harmonic crystal. It shows that in three di- 
mensions without introducing any pinning or disorder, 
harmonicity alone can give rise to a normal transport of 
heat in the crystal in the continuum limit. 
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